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PASSIVE  SONAR  DATA  PROCESSING,  PROPAGATION  MODELS 
AND  THE  DATA  CROSSCORRELATION  MATRIX  -  A  SURVEY 

by 

Even  B.  Lunde  and  Walter  M.  Zimmer 


ABSTRACT 


The  report  follows  one  possible  development  in  the  processing  of  data  from 
passive-sonar  arrays:  from  robust  processing  of  stationary  data  using 
simple  models,  to  refined  processing  of  data  acquired  under  non- stationary 
conditions,  using  more  complex  models.  For  such  a  development  to  take 
.place,  it  is  pointed  out  that  a  strong  interaction  between  researchers  in 
propagation  modelling  and  those  concerned  with  array  processing  seems 
Necessary.  As  the  level  of  ambition  related  to  information  extraction  from 
the  sound  field  is  increased,  some  of  the  most  common  assumptions  made  in 
the  processing  of  the  data,  such  as  stationarity,  plane  wavefront 
.propagation  model,  etc.,  have  to  be  questioned.  If  these  assumptions  do 
not  survive,  the  question  will  be  how  to  reformulate  them  and  how  to 
structure  new  models  in  such  a  way  that  they  can  be  incorporated  in  the 
data  processing  (information  extraction)  so  as  to  improve  sonar 
performance,  o 
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INTRODUCTION 


Let  us  start  out  with  some  definitions: 

X  :  Data  column  vector  with  N  elements,  where  the  elements  are  some 
representation  of  the  outputs  from  the  N  hydrophones  of  a  towed 
array,  being  part  of  a  passsive  sonar  system. 

R  :  E{X  X*}:  The  expected  cross-correlation  matrix  of  the  data  vector 
X.  Without  loss  of  generality,  it  is  assumed  that  E{X}=0.  E{  } 
is  the  expectation  operator,  i.e.  ensemble  averaging  operator.  *  is 
the  complex  conjugate  transpose  operator. 

In  a  surprisingly  large  part  of  all  passive  sonar  systems,  processing  on 
the  data  vector  X  involves,  explicitly  or  implicitly,^  the 
cross-correlation  matrix  R,  or  more  correctly,  its  estimate,  R. 

As  an  example,  a  weighted  summation  of  the  elements  of  X  is  often 
performed  (beamforming) 

y  =  ^w.*xi  =  W*X  , 

followed  by  an  averaging  of  the  squared  absolute  value,  the  power,  of  y: 

P  =  E{ |y |2}  E{yy*}  =  E{W*XX*W} 

=  W*E{XX*}W  =  W*RW 

Indeed,  nothing  but  the  cross-correlation  matrix  R  is  involved. 

The  two  main  reasons  for  the  extensive  concentration  on  the  R  matrix  are: 

The  probability  density  of  the  data  vector  X  is  often  assumed 
to  be  gaussian  (normal)  with  zero  mean  value.  Hence  the  density 
is  completely  described  by  the  cross-correlation  matrix  R  . 

Many  models,  criteria,  and  principles  used  in  the  data  processing 
are  such  that  only  first  and  second  order  statistics  are  asked 
for.  Typical  for  this  group  is  linear  least-mean- squares 
processing. 

Using  the  cross-correlation  matrix  in  array  processing  often  involves  two 
distinct  steps. 

Estimation  of  the  cross-correlation  matrix  itself. 

Using  the  matrix  to  extract  information  from  the  data  X  or 
derived  quantities,  like  R  itself. 

These  two  steps  are  usually  made  independently,  even  if  this  might  not  be 
optimal  [1]. 
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The  basic  assumption  in  estimating  R  is  stationarity  (or  quasi- 
stationarity) ,  such  that  instead  of  expectation  or  ensemble  averaging,  time 
averaging  can  be  used  as  a  substitute.  This  rather  crude  assumption  leads 
surprisingly  often  to  satisfactory  results.  The  reasons  must  be  that  in 
some  (few  ?)  cases,  the  process  in  question  is  indeed  stationary,  or  the 
employed  data  processing  is  robust  enough  to  tolerate  faulty  assumptions. 
In  the  present  report,  we  will  start  out  from  the  safe  position  of  robust 
processes  on  stationary  data  and  end  with  more  refined  modelling  where  lack 
of  stationarity  will  highly  influence  the  end  product  or  the  processing 
strategy  (as  it  happens,  this  is  also  the  chronological  order).  It  seems, 
also,  that  we  might  be  heading  for  a  future  abandonment  of  the  matrix, 

which  we  may  or  may  not  consider  to  be  disgraceful.  To  follow  this 

development  is  the  first  purpose  of  this  report. 

This  report  does  not  deal  explicitly  with  the  details  of  estimating  the 
matrix,  but  assumes  that  the  desired  estimate  is  available.  This  might  be 

of  the  matrix  itself,  R,  some  function  of  it,  such  as  the  inverse  R"1  [2], 
or  an  estimate  constrained  to  have  a  certain  structure,  typically  a 
Toeplitz  structure  [3].  The  explicit  estimate  itself,  R,  often  has  a  heavy 
influence  on  the  quality  of  the  ensuing  processing. 

Instead,  we  concentrate  on  the  task  of  extracting  information  from  the 
sonar  data  through  array  (spatial)  processing.  In  order  to  do  that,  we 

need  a  sound-propagation  model.  And  that  leads  directly  to  the  second 
purpose  of  this  report:  to  throw  some  light  on  the  relation  between 
propagation-model  research(ers)  and  array-processing  research(ers). 
According  to  [4]  the  two  disciplines  can  be  defined  respectively  as: 

-  To  improve  the  representation  of  the  sound  field  in  order  to 

predict  its  value  at  a  point  P  for  a  source  at  point  S  . 

-  To  refine  data  processing  in  order  to  discriminate  different 
sources  and  to  determine  their  individual  signals  together  with 
their  individual  position. 

It  goes  without  saying  that  a  strong  degree  of  interaction  between 
researchers  of  the  two  disciplines  is  highly  desirable.  This  should 
include  groups  acquiring  enough  interdisciplinary  knowledge  (and  interest) 
to  trigger  discussion  fruitful  to  both  disciplines.  To  improve  sonar 
performance  by  improving  data  processing  could  (and  should  ?)  be  the  main 
task  of  propagation  modelling.  Today  it  is  mainly  delegated  to  the 
(humble)  task  of  transmission- loss  predictions,  and  in  the  same  vein, 
evaluation  from  an  operational  point  of  view  of  a  new  sonar-system  concept. 
But  improved  propagation  modelling  is  not  often  part  of  the  new  concept. 
One  wonders  why  ! 

In  this  report  we  start,  as  an  example,  with  the  common  plane-wavefront 
propagation  model,  and  later  refine  it  to  the  so-called  normal-mode  model 
(5).  We  then  study  the  consequences  that  this  might  have,  both  on  the 
processing  itself,  and  on  the  amount  of  information  extracted  from  the 
data. 
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The  report  makes  no  claim  in  any  sense  to  be  complete  or  authoritative. 
Its  main  purpose  is  to  encourage  two  (independent)  lines  of  thought 
provoked  by  the  questions: 

-  How  healthy  are  some  of  the  assumptions  made  in  the  data 
processing?  Will  (and  should)  they  all  survive  as  more  refined 
models  are  applied  ? 

-  How  can  (and  should)  more  refined  sound-propagation  models  be 
formulated,  such  that  they  can  be  applied  to  the  processing  of 
sonar  data  so  as  to  extract  more  information  ? 
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1  ASSUMPTIONS 

1. 1  Frequency  bandwidth 

The  data  are  assumed  to  be  narrowband-f i 1 tered  with  bandwidth  B,  shifted 
down  to  zero  frequency  (quadrature  bandpass  filtering).  The  two  main 
reasons  for  this  are: 

A  delay  r  can  be  implemented  by  multiplication  with  a  phase 
rotating  factor,  exp{- j27tft} ,  as  long  as  the  delay/bandwidth 
product  is  much  smaller  than  unity,  t*B<<1  ,  where  f  is  the 
centre  frequency  of  the  band. 

If  a  stationary,  random  process  is  observed  for  a  time  that  is 
large  compared  with  its  correlation  time,  its  complex  fourier 
coefficients  are  uncorrelated. 

These  two  properties  often  simplify  implementation  -  such  as  implementing 
the  delay  necessary  in  beamforming,  where  the  maximum  delay  is  the  travel 
time  of  sound  across  the  array  -  and  allow  separate  treatment  of  the 
different  frequency  bands. 


1.  2  Receiver  system 

Our  receiver  system  is  a  towed,  horizontal  array  with  equidistant 
hydrophones  (Fig.  1).  The  distance  between  neighbouring  hydrophones  is 
Ax  and  the  total  number  is  N  . 

Both  spherical  and  Euclidean  coordinate  systems  are  introduced  (Fig.  1). 
The  first  hydrophone  (number  0)  is  the  reference  point,  except  for  the 
depth  coordinate,  which  is  zero  at  the  surface: 

x  =  R  cosp  y  =  R  sinp  cosiji  z  =  +  R  sinp  sin<J> 


(O.y.z) 


FIG .  1  RECEIVER  SYSTEM 
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2  PROPAGATION  MODELS 


1 


We  deal  with  two  models  only:  the  standard  plane-wavefront  propagation 
model,  and  the  normal -mode  model. 

2. 1  Plane-wavefront  model 

This  is  both  the  simplest  and  most  common  model.  It  is  the  tangential 
approximation  to  the  spherical  spreading  of  sound  from  a  point  source.  If 
this  tangential  plane  is  parallel  to  the  y-z  plane,  we  can  express  the 
sound  pressure  by: 


P(r,tn) 

-  ikr 

=  c  .  p  J 

n  e 

k  =  2n  A 

=  2nf/c  , 

where 

s  = 

stochastic  process  slowly  changing 

k 

wavenumber 

A 

wavelength 

f 

centre  frequency 

c  = 

sound  speed 

r  = 

distance  (displacement)  orthogonal  to  the  plane 
wavefront. 

The  plane-wavefront  approximation  is  assumed  valid  over  a  distance  d  in 
the  plane  satisfying,  A*r>d2,  the  so-called  far-field  condition,  where  r 
is  the  distance  from  the  source. 

Combining  this  model  with  the  receiver  array  of  Sect.  1.2  yields: 


«v 


P. 

g- jkx.cosp 

i  l 

i  n  i 

sn*D(k  cosp) 


where 

D  =  source  vector 

x.  =  position  of  hydrophone  i 

h  =  direction  of  incoming  plane  wavefront. 
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2. 2  Normal -mode  model 

The  normal-mode  model  |5]  is  a  natural  refinement  to  the  traditional 
plane-wave  model.  For  a  small  array  system,  the  modes  are 
indistinguishable,  while  they  will  be  more  and  more  noticeable  with 
increasing  array  size.  Using  the  simplest  possible  notation,  we  express 
the  sound  pressure  as: 


where 


M 


p(r,t  )  =  s  I  b_  e 
r  n  n  ,  m 
m=l 


-ik  r 
J  m 


k  =  horizontal  wavenumber 

m 

r  =  horizontal  displacement 

bffl  =  modal  weight,  being  a  function  of  source  depth, 
repeiver  depth,  and  range. 


Combining  this  with  our  receiver  array  yields: 


X(tn)  =  s 


M  (  -jk  x. 
n  I  b  l  e  m  1 
nm=l  m  ( 


cosp)  M 

=  s  I  bm  D(k  cosp)  , 


/  j  t— 

;  nm=l 


m  m 


sni,b»  D<k  cosbm^  =  D(V  • 


m=l 


where 


cos^m  =  T  cosp  . 


This  expression  says  that  the  normal-mode  model  can  be  interpreted  as  a 
weighted  combination  of  plane  waves,  each  with  slightly  different 
direction,  dependent  on  the  mode's  horizontal  wavenumber,  kffl  .  It  is  worth 

noticing  that  km  is  a  function  only  of  the  local  channel  (medium). 
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3  PROCESSING  WITH  PLANE-WAVEFRONT  MODEL 

In  this  chapter  we  study  the  connection  between  the  plane-wavefront  model 
and  the  cross-correlation  matrix.  It  is  natural  to  introduce  the  concept 
of  sound-field  density.  Given  the  matrix,  what  information  is  available 
about  the  sound-field  density,  and  can  it  be  estimated  uniquely  ? 
Conventional,  adaptive,  and  maximum-entropy  beamforming  give  different 
answers.  We  will  see  why. 


3. 1  Cross-correlation  matrix  and  the  sound  field 

According  to  Sect.  2.1,  the  array  datum  from  one  point  source  is: 

X  =  sn*D(k  cosp )  .  (Eq.  1) 

The  contribution  to  the  cross-correlation  matrix  would  be: 


R  =  E{XX*}  =  E{sn-sn*  D(p)-D*(p)}  (Eq.  2) 

=  E{sn-sn*}  D(p) •  D*(p)  =  a2(p)-D(p)-D*(p)  , 

where  E{  }  is  the  expectation  operator  and  s^  is  assumed  to  have  zero 

mean.  Equation  2  can  be  generalized  to  include  independent  contributions 
from  al 1  directions: 


o2(P)  -*  P(p,4i)  sinp  dp  di|i  (Eq.  3) 

R  =  J  dp  J  dip  P(p,40  sinp  D(p)  D*(p)  ,  (Eq.  4) 

P 


where  P(p,t|i)  is  the  sound-field  density.  Obviously,  as  D(p)  is 

independent  of  <|j  ,  the  horizontal  line  array  integrates  over  tji  : 


where 


R  =  X  dp  si np  P(p)  D(p)  D*(p) 

P 


n 

P(P)  P(p,4>)  ^ 

0 


(Eq.  5) 


(Eq.  6) 
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chosen  so  that  the  model  matrix  R  fits  as  closely  as  possible  to  the 
estimated  matrix  R  (from  data).  As  R  has  a  Toeplitz  structure,  R  can 
be  estimated  while  being  constrained  to  such  a  structure  [3].  In  that 
case  R  provides  2N-1  real  data  values: 


r^  ,  I k I  g  N ,  wi th  r^  =  r* 


Hence,  to  determine  all  the  parameters  requires  that  K  +  T  ^  N 

One  approach  to  this  estimation  problem  is  to  attempt  to  estimate  the 
continuous  part  first.  As,  by  definition,  the  discrete  part  should  produce 
a  matrix  R^  of  rank  T  %  N-K  ,  the  parameters  a^,  I k |  <  K  ,  should  be 

chosen  such  that  the  K  smallest  eigenvalues  of  R  ~  Rc  are  as  close  to 
zero  as  possible.  The  estimation  criterion  in  [9]  is  in  agreement  with 
thi s  pri nci pi e: 


N-l 

min  J  =  1  K2 

B  n=N-K  " 


(Eq.  47) 


K-l  N-l 

R  =  R-R  (B)  =  R  -  I  b  I  =  2  A  E  E*  (Eq.  48) 

a  c  k=l-K  K  n=0 

where 

B  =  [bQ,  br  -  b^]1  (Eq.  49) 

An=  Eigenvalues  of  R  -  Rc(B) 

En=  Eigenvectors  of  R  -  Rc(B) 

Details  about  the  algorithm  are  given  in  App.  A. 

The  next  step  is  to  estimate  a2,  and  v.,  (i=0,  ...,  T-l)  from  R  ,  which 

11  3 

is  now  identical  to  the  same  problem  in  the  Pisarenko  method;  again  the 
goniometer  method  is  one  candidate  [8]. 

If  we  apply  this  method  to  our  old  example  of  one  point  source  and  white 
noise,  we  of  course  obtain  the  correct  answer  (our  example  is  chosen  to  fit 
the  model  !  ) 

3. 2  Sensitivity  to  Source-Vector  Model 

All  the  methods  outlined  in  the  previous  sections  assumed  a  plane-wave 
propagation  model,  leading  to  the  classical  steering  vector  as  a  model  of 
the  source  vector. 
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Comparison  with  Eq.  36  yields  T  =  N-l.  At  this  stage,  several  methods  can 
be  used  to  estimate  a?  and  v.#  i  =  1,N-1.  Among  them  is  the  so-called 

goniometer  method  |8J. 


In  order  to  allow  a  more  general  continuous  term,  the  authors  have 
attempted  to  generalize  the  Pisarenko  method  [9].  One  possibility  is  to 
view  Eq.  37  as  the  zero  term  in  a  fourier  series.  A  generalized  version  is 
then  obviously: 


pc  ( v  ) 


such  that 


b 


k ,  c 


K-l 


k  1  -  K 


- j2nkv 


-  a,  ,  I  k  I  <  K 
k 

=  0  .  Ik!  ^  K 


The  matrix  is: 


(Eq.  42) 


K-l 

R  =  2  a  I, 

c  k=l-K  k  k 


where  is  the  zero  matrix,  except  that  diagonal  number  k  has  elements 

of  value  one.  For  K  =  1  this  corresponds  to  Pisarenko  method. 

To  summarize,  the  structure  is: 


P(v) 


T-l 

2 

i=o 


a?  6( v-v . )  + 
7  7 


K-l 

2 

k=l-K 


\  e 


- j2nkv 


(Eq.  43) 


T-l 

2 

i  -o 


i  ej2"kv<  *  \ 


b 


k 


=  < 


T- 1  .„  . 

I  o?  eJ27lkvi 
i=o 


1  k  1  <  K 


I  k  |  £  K 


(Eq.  44) 


(Eq.  45) 


R 


T-l  K-l 

1  a2.  D( v  . )  D*( v . )  +  1  a.  I. 

i=o  k=l-K 


(Eq.  46) 


where  the  first  term  in  all  these  expressions  is  the  discrete  one  (signal 
part)  and  the  second  term  is  the  continuous  one  (noise  part).  Altogether 
there  are  2K-1  +  2T  unknown  real  parameters  (a^  =  a*  _^).  They  have  to  be 
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It  is  quite  important  in  the  following  that  the  rank  of  is  given  by: 

rank  (Rd)  =  min(T.N)  ,  (Eq.  36) 

where  N  is  the  dimension  of  R^  and  the  number  of  hydrophones.  The 

orthogonal-decomposition  method  will  work  only  if  this  rank  is  less  than  N; 
in  other  words,  if  R^  is  a  singular  matrix.  This  means  that  only  T  <  N 

discrete  lines  can  be  allowed  in  the  model  of  P(v). 

Turning  attention  to  the  continuous  term  P  (v),  it  is  not  so  obvious  or 
straight-forward  how  to  model  Pc  or  Rcc  .  Mermoz  [7]  indicates  that 

general  plausible  models  of  noise  turbulence  could  be  used  to  formulate  Rc 
as  a  function  of  several  free  parameters  A  (A  is  a  parameter  vector): 

R  =  R  (A)  . 
c  cv 

and  that  A  could  be  calibrated  (chosen)  in  such  a  way  that  the  matrix 
R^  =  R  -  Rc  has  minimum  rank.  If  the  rank  is  larger  than  the  number  of 

free  parameters  in  Rc  ,  it  indicates  a  good  noise  model.  Bienvenu  and 

Kopp  |8]  work  along  these  lines. 

Another  favourite  is  the  Pisarenko  model,  which  has  only  one  parameter: 

Pc(v)  =  cr£  .  (Eq.  37) 

The  motivation  is  that  Pc(v)  contains  uncorrelated  noise  in  the  array 
data.  Note  that  in  this  case: 

Rc  =  ^  I  ,  (Eq.  38) 

T-l 

R  =  Rri  +  R  =  I  a?  D(v.)  D*(v.)  +  o2  I  ,  (Eq.  39) 

L  i=o 

or 

Rd  =  R  -  Rc  =  jj  <\  -  °P  En  E„*  •  <E“-  40> 


where  and  En  are  the  eigenvalues  and  eigenvectors  of  R  . 

In  order  to  give  Rd  a  rank  less  than  N  ,  is  chosen  to  be  equal  to 
the  smallest  eigenvalue  of  R  ,  An-  resulting  in: 


N-l 

1  (A  -  o2)  E  E  * 
n=l  n  c  n  n 


(Eq.  41) 


rank  (Rd)  =  N-l 
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It  is  worth  mentioning  that  the  constraints  (Eq.  30)  can  be  disguised,  or 
made  somewhat  different,  without  changing  the  principle.  The  so-called 
WB2  algorithm  uses  the  following  constraints  [6]: 


P(uk)  =  D*(uk)  R  D(uk)  = 


-u 


dv  P(v)  |D*(uk)-D(v)j  2 


(Eq.  31) 


This  means 
"beampattern 
of  Eq.  31. 


that  the  WB2  algorithm  has  replaced  exp{-j2nkv}  by  the 
"  |D*(uk)  D(v)|  2  as  the  weighting  function  in  the  functional 


3.1.5  Continuous  plus  discrete  decomposition 
(Orthogonal  decomposition  method) 

Because  of  the  limited  size  of  the  array,  all  the  previous  methods  have 
some  problems  with  estimating  the  discrete  point  source  in  the  sound-field 
density.  One  way  to  improve  that  is  to  allow  for  a  finite  number  of 
discrete  lines  in  the  estimation  of  the  density,  such  that  P(v)  consists 
of  a  discrete  and  a  continuous  term: 

P(v)  =  Pd(v)  +  Pc(v)  .  (Eq.  32) 


Keeping  the  plane-wave  model,  the  discrete  term  is  easy  to  structure 
T-l 

P.(v)  -  2  •  Mv-v.)  ,  (Eq.  33) 

u  i=0 


where  T  is  the  number  of  sources. 


Ihe  contribution  to  the  fourier-ser ies  coefficient  is: 


b 


k,d 


/  Pd(v)  eJ2nkv  dv 


T-l 


j2okv  ■ 
e 


(Lq.  34) 


and  to  the  cross-correlation  matrix  is 


h  T-l 

Rh  =  J  P.(v)  D(v)  D*(v)  dv  =  /  <i-  l)(v  )  ■  D*  ( v  _  (tq.  3b) 

°  -k  0  (0 


t 
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Hence  this  method  seems  to  give  good  estimates  both  for  very  low  and  for 
very  high  c|2/o2n  ratios,  and  uses  the  available  information  consistently. 

Therefore 


b  a  Ikl*  N  , 


while  it  is  somehow  more  imprecise  in-between.  However  it  is  worthwhile 
noticing  from  Eq.  28  that  to  obtain  the  peak,  the  "propagation  model"  D(v) 
has  to  match  perfectly  the  "real  propagation",  0(o)  ,  when  v  =  0 

Especially  is  this  the  case  for  high  values  of  a2/o2. 


3.1.4  Maximum-entropy  method 

So  far,  except  for  the  truncation  method,  none  of  the  methods  have 
explicitly  bothered  about  using  the  available  ^information  in  a  consistent 
way,  in  the  sense  of  explicitly  requiring  "6^  =  b^  for  Ikl  S  N.  The 

maxi  mum- entropy  method  does  that.  The  other  fourier-series  coefficients 
are  then  determined  according  to  the  maximum-entropy  principle  [3]: 


max  J  =  J"  dv  In  P(v)  ,  (Eq.  29) 

p(») 

constrained  by: 

H  .  o 

rk  =  ;  P(v)  e"j2nkv  dv  for  Ikl  <  N  .  (Eq.  30) 

Also  in  this  case  an  analytical  expression  has  not  been  obtained  for  b^ 

when  lk|>N.  Computer  plots  (Fig.  5)  show  for  different  S/N  ratios 
^a  =  0  dB,  b  =  -10  dB,  c  =  -20  dB)  the  expected  behaviour  that 
bk  =  bk  =  r^  for  I k I <  N,  and  then  a  smooth  decrease  in  bk  approaching 

zero  for  large  values  of  k  ,  in  agreement  with  the  maximum-entropy 
principle,  (or  "minimum  a-priori  assumption  principle"). 

/•v 

The  P(v)  plots  show  excellent  "resolution",  as  the  peak  is  even  narrower 
than  the  maximum  likelihood  peak.  However,  the  warning  about  sensitivity 
is  even  more  appropriate  in  this  case. 
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The  result  is: 


.  1 


_  i 


P(v)  -  |D*(v)  R  0( v ) ]  =  J  du  P(u) 


D*(v)R~  D(u)  12 

D*(v)R~  0( v) 


In  our  example: 


N  a2  +  a2 
s  n 


P(v) 


N  as  ..  I  D*(v)  D(o) 

"a7"  N 

n  L  1 


1  + 


(Eq.  28) 


The  corresponding  is  not  found  analytically,  but  is  shown  in  Fig.  4 
together  with  p(v)  as  a  function  of  S/N  =  cr2/o2  for  N  =  32  . 


The  figure  shows  that  the  method  assigns  a  value  ^  0  for  |k|  Z  N,  in 

such  a  way  that  there  is  a  smooth,  rounded  transition  near  |k|  =  N, 
depending  on  o£/o*. 


This  results  both  in  low  or  almost  no  oscillations  (no  sidelobes),  and  a 
narrow  peak  at  the  position  of  the  source.  Hence  it  is  a  big  improvement 
over  conventional  beamforming,  in  the  sense  that  it  is  a  pre-whitened 
matched  filter.  It  takes  into  account  the  actual  sound-field  distribution, 
to  the  extent  that  the  matrix  R  can  11  image"  it.  The  disadvantage  is  of 
course  that  the  available  information  is  still  not  used  in  a  consistent 
way,  because  i  b^  for  k  <  N. 

To  gain  some  insight,  we  consider  some  limiting  cases,  a2  -*  0 
for  fixed  value  of  cr2  The  first  case  is  trivial,  because 
al  I  values  of  k  : 


k  =  0 
k  *  0 


As  o2s>>  a2  ,  the  peak  in  the  P(v)  plots  becomes  both  higher  and 
more  narrow  (at  -3  dB  below  peak  level  as  an  example)  with  increasing  a2s  . 

The  plots  show  that  the  transition  from  to  b^»0  happens 

both  much  slower  and  at  increasingly  higher  k  levels. 


and  a2  ->  ® 
s 

bk=bk ,or 
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As  seen,  this  method  is  using  the  available  information  in  a  consistent 
way,  1)^  =  for  |k|<  N,  but  the  effect  of  truncation  creates  heavy 

oscillations  in  the  estimate  of  P(v).  It  might  even  result  in  some 

negative  estimates  of  P(v)  for  some  values  of  v  ,  depending  on  the 

"signal-to-noise"  value,  a2  /  a2  .  Hence  the  method  has  certain 

s  n  , 

I /N  J 

disadvantages.  Figure  2  shows  a  plot  of  |b^|  and  P(v)  for  different 
signal-to-noise  ratios. 


3.1.2  The  conventional  beamformer  (the  matched  filter) 


The  expected  normalized  power  output  of  conventional  beamforming  is  given 
by: 


P  (v)  D*(v)  R  D(v) 


=  a2s  •  jjj  |D*(v)  •  D(o)  | 2  +  a2 

=  a2  .  I  fsin(7i  N  v)~|2  +  2 

as  N  sin(n  v)  J  an 


b  k  =  (N  -  Ik!) 
=  0 


1  k  I  <  N 
I  k  |  >  N 


(Eg.  24) 


(Eq.  25) 


As  seen,  this  method  not  only  truncates  the  fourier  series  but  also  weights 
the  known  values  with  a  triangular  window,  the  so-called  Bartlett  window. 
This  has  the  advantage  of  avoiding  high  oscillating  levels  (sidelobe 
levels),  but  at  the  same  time  broadens  the  peak  of  our  point  source  (giving 
a  broad  main  lobe).  Another  disadvantage  of  the  weighting  is  an 
inconsistent  use  of  available  information.  As  a  last  comment,  this  method 
always  yields  non-negative  estimates  of  sound-field  density.  Plots  of  |b.| 
and  P(v)  are  given  in  Fig.  3  for  different  signal-to-noise  ratios. 


3.1.3  Maximum-likelihood  beamforming 

The  name  of  this  method  is  misleading,  as  it  is  in  reality  a  minimum 
expected-power  beamforming  method,  constrained  by  constant  gain  in  the 
beam-steering  direction: 


P(v)  =  min  W*(v)  R  W(v)  =  min  f  du  P(u)|w*(v)  D(u)| 2  ,  (Eq.  26) 
W(v)  W(v)-*s 

constrained  by: 

W*(v)  D(v)  =  constant  =  N  (Eq.  27) 
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From  Eqs.  18  and  19,  we  can  say: 

-  Px(v)  is  completely  determined  by  R,  and  vice  versa. 

-  P2(0  is  completely  independent  of  R,  and  vice  versa. 

-  As  P(v)  =  Pj(v)  +  P2(v),  the  sound  field  density  can  never  be 
uniquely  determined  from  the  cross-correlation  matrix  R  alone. 

This  means  that  to  obtain  a  unique  estimation  of  P(v)  ,  either  some 
additional  information  is  needed  or  some  estimation  criteria,  constraints, 
or  principles  have  to  be  applied.  This  is,  of  course,  the  reason  that  we 

have  such  a  jungle  of  estimation  methods  and  proposals,  more  or  less 

consistent  with  the  real  facts. 

In  the  following  sections  we  will  study  four  methods  for  estimating  P(v), 
considering  the  example  of  a  point  source  in  white  noise: 

P(v)  =  6(v)  +  ,  (Eq.  20) 

b,  =  J  P(v)  e^n^V  dv  =  a2  +  ,  k  =  0  (Eq.  21) 

K  i  s  n 

=  a2  .  k  t  0 


3.1.1  Truncated  fourier  series  method 

The  truncated  method  uses  the  available  information  about  the  fourier 
series  of  P(v)  and  neglects  the  rest: 

b.  =  b,  -  o2  +  a2  ,  k  =  0 

k  k  s  n  ’ 

k  *  0  ,  | k| <  N  (Eq.  22) 

|k|§  N 


The  corresponding  estimate  of  P(v)  is: 

P(V)  =  Pl(v)  =  V  b.  e'j2nkv 
k=l-N  K 


=  „2  .  slnJn.i2IM.J_ vl  +  2 
°s  sin(n  v)  °n 


(Eq.  23) 
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Substituting  Eg.  11  into  Eq.  10  yields: 


rk  =  J  dv  l  bg  e''i2n£v  ej2nkv  =  bk  .  |k|<  N 

,  £=-oo 

1 2 

This  shows  that  the  elements  of  diagonal  k  is  equal  to  the  k- 
coefficient  of  the  sound-field  density. 

To  see  the  full  consequence  of  this,  we  split  P(v)  into  two  terms 
P(v)  =  Pi(v)  +  P2(v)  , 

in  which 


(>,<»)=  "i1  b.  e-J2"'1'  , 

£=1-N  * 


P2(v)  =  Ib.e 
| £| SN  * 


-j2n£v 


From  this  follows: 

h 

Rj  =  J  dv  Pt(v)  0(v)  D*(v)  =  R 

~h 

h 

R2  =  J  dv  P2(v)  D(v)  D*(v)  =  0 

or  in  words: 

-  the  matrix  R  is  completely  determined  by  Px(v). 

-  the  matrix  R  is  completely  independent  of  P2(v). 


Also,  by  combining  Eqs.  13  and  15: 


P^v) 


N-l 

I 

£=1-N 


e*j2n£v 


11 


(Eq.  13) 

th  fourier 

(Eq.  14) 

(Eq.  15) 

(Eq.  16) 

(Eq.  17) 

(Eq.  18) 

(Eq.  19) 
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For  our  linear  equi-spaced  array  it  is  convenient  to  change  the  variable 

2n*v  =  -k*a*cosp  ,  C Eq .  7) 

which  yields: 
a/X 

R  =  ;  dv  P(v)  D(v)  D*(v)  .  (Eq.  8) 

-a/X 

We  also  include  the  virtual  field,  thereby  obtaining: 

h 

R  =  J  dv  P(v)  D(v)  D*(v)  .  (Eq.  9) 

R  is  always  a  positive,  semi-def ini te,  Hermitian  matrix 
R*  =  R  (Hermitian) 

h 

y*Ry  =  J  dv  P(v) I y*D( v)| 2  £  0  (positive  semi-definite)  . 

-h 

For  a  linear,  equidistant  array,  R  defined  by  Eq.  9  is  also  Toeplitz: 


m,n 


h 

=  /  dv  P(v) 


e-j27tmv  e+j27tnv 


r 

m-n 


kj  <N 


(Eq.  10) 


where  k  =  0  indicates  an  element  on  the  main  diagonal;  k  =  N-l  is  the 
element  in  the  lower  left  corner  and  k  =  1-N  is  the  one  in  the  upper  right 
corner  of  the  R  matrix. 

We  now  turn  to  a  quite  important  problem,  judging  by  the  numbers  of  papers 
published  in  the  open  literature:  how  well  can  we  estimate  the  sound-field 
density  P(v)  from  a  knowledge  of  R  ?  This  is  exactly  the  same  problem 
as  estimating  the  power  spectrum  P(f)  from  an  estimate  of  the 
autocorrelation. 

To  see  this,  we  write  P(v)  as  a  fourier  series: 


oo 

P(V)  =  I  b. 

£=-t»  * 


e-j2n£v 


> 


where,  since  P(v)  is  a  real  quantity, 


(Eq.  11) 


(Eq.  12) 
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We  know  that  this  model  is  not  correct  in  the  real  stratified  ocean.  It 
seems  to  be  quite  a  good  model  for  modestly  sized  arrays,  but  as  soon  as 
the  ambition  for  higher  resolution  leads  to  both  larger  arrays  and  higher 
resolution  methods,  one  might  be  heading  for  problems. 

Improved  propagation  models  have  clearly  shown  that  the  horizontally 
layered  propagation  channel  consisting  of  the  surface,  ocean,  bottom,  and 
sub-bottoms,  gives  rise  to  several  modes,  each  having  different  horizontal 
wavenumbers.  As  soon  as  this  difference  is  sensed  by  the  combination  of 
sufficiently  larger  array  and  sufficiently  high-resolution  method,  the 
mismatch  between  reality  and  the  plane-wavefront  model  might  cause  serious 
errors  in  parameter  estimation. 

In  the  next  section  we  therefore  study  what  can  be  done  with  the  normal¬ 
mode  model  in  order  to  improve  the  situation,  and  we  address  eventual 
problems  in  using  it. 
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4  NORMAL-MODE  PROPAGATION  MODEL 
4. 1  Principles 

In  this  section  we  substitute  the  plane-wavefront  source  vector  D(p), 
with  the  compound  normal -mode  source  vector: 

M 

F(p)  =  I  a  D(p  )  ,  (Eq.  50) 

m=l 


where  8  is  defined  by: 
m 

k  cos  0  =  k  cos  p  (Eq.  51) 

m  m 

One  might  think  that  it  is  simple  to  use  F(p)  instead  of  D(p)  in  all 
the  traditional  methods.  But  there  is  one  important  difference.  While 
D(p)  depends  on  only  one  environmental  parameter,  which  is  the  local  sound 
speed  around  the  array,  F(p)  has  a  much  more  complicated  dependence  on 
the  channel  parameters.  In  fact,  in  order  to  know  the  mode  weight  a^,  or 
more  generally,  the  mode  vector: 

A  -  [  aj  . . . ,  a^  j  , 

one  has  to  know  not  only  the  local  normal-mode  characterization  of  the 

channel,  but  the  channel  parameters  all  the  way  from  the  source  to  the 

receiver.  It  is  a  highly  unrealistic  assumption  that  the  mapping  of  the 

ocean  and  the  bottom  parameters  will  ever  be  so  detailed  as  to  allow  an 
a  priori  calculation  of  A  with  the  necessary  accuracy.  Imagine  only  the 
problem  of  obtaining  a  satisfactory  value  for  the  relative  phase  of  the  A 
elements.  To  do  so,  we  really  need  to  know  the  range  to  the  source  quite 

well.  As  to  the  value  of  the  relative  amplitude,  this  depends  not  only  on 

the  relative  attenuation  of  the  modes,  but  also  on  the  depth  of  the  source. 
And  this  even  neglects  the  effect  of  any  randomness  in  the  channel,  or  any 
movement  of  the  source  or  receiver. 

The  above  arguments  have  convinced  the  authors  that  the  model  vector  A  can 
never  be  calculated  satisfactorily  from  scratch  (channel  parameters);  hence 
it  is  not  possible  to  obtain  an  a  priori  source  vector  model  F(p).  So  far 
the  only  information  we  have  is  a  structure  for  F(p),  given  by  Eq.  50. 

Looking  to  the  modal  source  vectors  D(0m),  they  depend  on  the  horizontal 

wavenumbers  k^  .  But  as  the  source  vectors  really  express  relative  phase 

rotation  (or  phase  delay),  it  is  the  local  horizontal  wavenumbers  that  are 
required.  Hence  it  does  not  really  matter  if  the  channel  is  range- 
dependent  between  source  and  receiver,  as  long  as  we  know  the  local  channel 
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parameters  quite  well.  We  will  make  the  following  assumption: 

-  The  local  envi ronmental  parameters  are  known  so  well  that  the 
horizontal  wavenumbers  can  be  calculated  with  sufficient 
accuracy. 

What  we  mean  with  "sufficient  accuracy "  depends  on  the  array  size  and 
processing  method. 

At  this  point  we  have  obtained  a  structuring  of  F(p),  according  to 
Eq.  50,  in  which  all  the  source  vectors,  D(pm),  except  the  bearing,  are 

known  but  where  the  mode  weights,  a  ,  are  unknown  and  have  to  be 

estimated.  On  one  hand  this  gives  many  more  parameters  to  estimate  than 
for  the  plane-wavefront  model.  On  the  other  hand,  it  results  in  a 
non-Toeplitz  matrix  containing  much  more  information. 

It  seems  obvious  that  methods  requiring  a  priori  knowledge  of  the  source 
vector  are  useless  in  this  case.  Hence,  we  concentrate  on  the  orthogonal 
decomposition  method  in  the  next  section. 


4.  2  Orthogonal  Decomposition  Method 

It  still  seems  convenient  to  split  the  cross-correlation  matrix  into  two 
terms,  one  caused  by  discrete  point  sources  and  one  by  continuous 
(distributed)  ones,  such  as  the  sea  surface.  The  discrete  term  is: 


Rd  = 


T-l 

I 

i=0 


a2.  F .  (B . )  F.*(p.)  , 

l  i  Vhi'  ’ 


T-l 

=  1  a2  G(p.)-A.-A.*-G*(B.) 

i=0  iii  i 


(Eq.  52) 


where  we  have  introduced: 


M-l 


F  .  (p  . )  = 

=  I  a  . 

m,  l 

D(p  .) 
Hm,  i  9 

-  GO,) 

*Ai 

(Eq. 

53) 

n=o 

G. 

l 

=  [D  .... 
0,1 

D 

m,  i 

W 

] 

(Eq. 

54) 

Ai 

=  [a  .... 

1  o,i 

a 

m ,  i 

aM-r' 

)T 

(Eq. 

55) 

in  which  A. , 

l  * 

a?  and  p. 
i  *i 

are  the  unknown  quantities 

to  be  determined. 

Notice  that  the 

rank  of  the 

matrix 

is  T  . 

every  source 

contributing 

wi  th 

rank  1. 


A  similar  model  of  the  continuous  part  would  be  forbiddingly  complex. 
Luckily,  in  this  case  the  intermodal  terms  have  a  tendency  to  cancel  each 
other  [10).  Therefore  the  structure  again  simplifies  to  one  that  is 
identical  to  the  plane-wavefront  case: 


Rc  =  J  dv  P  (v)  D(v)  D*(  v) 


K-l 

*  bk 
k=l-K  ’ 


I 


k 


(Eq.  56) 


The  total  model  is  then: 


T-l 

R  =  I  erf  G(p.)  Ai  A.*  G*(p.)  + 


K-l 

Z  bk 
k=l-K  *• 


(Eq.  57) 


Hence  we  can  again  use  the  same  estimation  criterion,  Eqs.  47  and  48,  and 
the  same  algorithm  as  in  the  plane-wavefront  case  of  Sect.  3.1.5.  The  rank 
of  R^  does  not  change  with  the  change  of  model  as  long  as  the  modal  source 

vectors  add  coherently. 

The  problem  is  now  to  estimate  the  unknown  parameters  of  R^  by  matching  R^ 

and  R  .,  both  assumed  to  have  reduced  rank  T<N: 
d 

One  procedure  to  determine  those  parameters  is  given  in  App.  8. 
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5  NORMAL-MODE  PROPAGATION  MODEL  WITH  DYNAMIC  SOURCE/RECEIVER  GEOMETRY 

As  long  as  the  source  and  the  receiver  do  not  move,  we  have  assumed  that 
there  are  no  problems  in  obtaining  a  good  estimate  of  the  cross-correlation 
matrix,  using  time  averaging.  The  introduction  of  dynamic  source/receiver 
geometry,  however,  introduces  an  essential  difference  between  the  classical 
plane-wave  source  vector  (delay  vector)  and  the  normal-mode  source  vector. 

In  the  plane-wave  case,  the  source  vector  normally  changes  quite  slowly 
when  the  source  is  in  the  far  field.  Hence  it  might  still  be  possible  to 
obtain  quite  a  good  estimate  of  the  matrix,  given  that  the  receiver  system 
does  not  change  course. 

It  is  quite  a  different  situation  with  the  normal-mode  source  vector.  As 
mentioned  earlier,  this  consists  of  a  weighted  sum  of  source  vectors  for 
the  individual  modes.  These  individual  source  vectors  behave  exactly  like 
the  plane-wave  source  vector.  The  relative  phases  of  the  weights,  however, 
rotate  with  time,  causing  the  individual  vectors  to  combine  into  a 
resulting  source  vector  that  changes  quite  drastically  with  time.  The 
deterministic  relative  phase  rotations  in  the  weights  have  two  effects: 

The  modes  have  different  horizontal  wave  vectors,  such  that  a 
displacement  in  space  results  in  relative  phase  rotation. 

The  modes  have  different  doppler  shifts,  such  that  displacement 
in  time  results  in  relative  phase  rotation. 

The  end  result  is  that  the  stationary  condition  for  replacing  ensemble¬ 
averaging  with  time-averaging  is  not  well  satisfied.  There  seem  to  be 
several  plausible  ways  to  resolve  this  dilemma.  Three  of  them  are  outlined 
in  the  next  sections. 


5 .  1  Update  of  Matrix  with  Stationary  Model 

We  have  seen  that  the  normal-mode  source  vector  changes  with  time  in  a 
dynamic  source/receiver  situation.  Let  us  see  what  happens  to  the 
contribution  to  the  vector  from  such  a  source  if  we  s  t i 1 1  prefer 
straightforward  t ime-averagi ng. 

To  repeat,  the  data  vector  X  can  be  written: 


R  -  ss* 


2 

m 


la  a 
m  n 
n 


m  n 
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s  is  independent  of  the  weights  a  and  D  is  assumed  not  to  change  over 
the  averaging  time: 


Because  of  the  relative  phase  between  the  weights,  a  a-*  «  0  for 
m  ^  n,  and  m  n 


This  expression  corresponds  to  M  independent  point  sources  using  a 
plane-wave  model . 

If  the  array  is  large  enough  or  has  enough  vertical  extension  to  resolve 
the  mode,  this  matrix  model  might  be  useful.  What  is  lost  is  information 
about  the  weight  phase,  which  is  sometimes  of  interest.  Nevertheless  the 
approach  is  simple  and  might  be  given  further  thought.  With  several 
sources  present,  however,  it  might  be  difficult  to  identify  the 
contributions  from  the  different  sources. 


Update  of  Matrix  with  Non-stationary  Model 


A  common  way  to  handle  a  non-stationary  situation  is  to  use  a  state/space 
approach,  which  in  the  linear,  gaussian  case  results  in  the  well-known 
Kalman-Bucy  filter.  Even  in  many  non-linear  and  non-gaussian  cases  it  has 
been  used  with  some  success.  What  are  required  in  this  approach  are  a 
state  model  and  a  measurement  model.  Very  little  work  has  been  reported  on 
the  use  of  such  a  method  to  estimate  (track)  the  cross-correlation  matrix 
itself  [11],  but  a  straightforward  Kalman-Bucy  filter  approach  does  not 
seem  feasible.  Some  work  is  being  done  [12]  on  the  use  of  an  input/output 
formulation  (the  opposite  of  the  state/space  formulation),  the  so-called 
scattering  framework  for  estimation,  to  handle  non-stationarities  but 
results  have  not  yet  been  published. 


5.  3  Abandonment  of  Matrix 

When  dealing  with  zero-mean,  gaussian,  stationary,  stochastic  processes  in 
linear  systems,  the  cross-correlation  matrix  yields  all  the  available 
information. 

Sections  5.1  and  5.2  indicate  that  maintaining  an  approach  based  on  the 
cross-correlation  matrix  might  lead  to  difficulties  in  dynamic  situations. 
This  is  even  more  so  if  non-linearities  and  non-gaussian  processes  enter 
the  picture.  Therefore  an  increasing  amount  of  attention  is  being  given  to 
more  general  approaches  capable  of  handling  quite  complex  situations  and 
many  books  and  conferences  are  partly  or  completely  occupied  with  these 
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problems  [13-16].  As  a  result  of  its  success  in  the  linear  gaussian  case, 
the  Kalman-Bucy  filter,  most  of  the  work  is  based  on  the  state/space 
approach. 

Work  in  this  direction  initially  led  to  large  problems  dealing  with 
non-linear,  stochastic,  partial  differential  equations  (PDE).  This 
situation  has  improved  drastically  in  recent  years,  going  through  an 
intermediate  phase  of  linear  stochastic  PDE  to  the  latest  progress  in  which 
a  linear  deterministic  PDE  is  parametrized  by  the  stochastic  process 
represented  by  the  measurements  [17]. 

SACLANT  ASW  Research  Centre  has  also  done  some  work  in  this  general 
direction  for  the  special  case  of  target-  motion  analysis  [18-20]. 

It  might  be  that  to  gain  maximum  benefit  from  improved  propagation  models, 
of  which  the  normal  mode  is  an  important  example,  one  has  to  turn  to  such 
general  methods. 
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CONCLUSION 


It  is  quite  common  in  passive  sonar  data  processing  to: 

-  Apply  a  plane  wavefront  propagation  model. 

“  Assume  stationary  processes. 

-  Let  the  cross-correlation  matrix  of  the  array  data  play  an  important 
part  in  the  processing. 

For  a  horizontal  array  of  modest  size  the  three  points  above  are  quite 
reasonable.  Shallow-water  sea  trials  with  towed  arrays  have  demonstrated 
that  the  wavefront  model  seems  almost  perfect  for  such  arrays.  The 
combination  of  modest  size  and  a  plane-wavefront  model  makes  the  processing 
quite  robust  with  respect  to  the  stationary  assumption.  Hence  the 
cross-correlation  matrix  can  be  estimated  quite  well  without  major 
problems. 

The  main  problem  with  a  small  array  is  its  limited  resolution.  The 
sound-field  density  can  be  estimated  only  with  fundamental  uncertainties, 
in  spite  of  high-resolution  methods.  The  necessary  information  is  not  in 
the  data. 

Increasing  the  array  size  leads  to  a  chain  reaction  of  consequences.  The 
first  consequence  might  be  that  a  pi ane-wavef ront  model  is  inconsistent 
with  the  data.  A  more  refined  propagation  model  is  then  necessary.  At 
least  for  low-frequency,  shallow-water  conditions,  a  normal-mode 
propagation  model  seems  a  natural  choice,  and  was  attempted  for  use  in  this 
report. 

This  leads  directly  to  the  second  consequence,  that  better  knowledge  of  the 
environmental  parameters  is  necessary.  The  plane-wavefront  model  required 
only  the  knowledge  of  the  sound  speed  around  the  array.  As  a  minimum 
(light  modelling)  the  normal-mode  model  requires  the  environmental 
knowledge  needed  to  calculate  the  modal  horizontal  wavenumbers  [7].  In 
fact  we  may  never  know  the  channel  well  enough  to  apply  some  traditional 
processing  schemes.  For  instance,  with  conventional  beamforming  in  which 
the  steering  vector  of  the  plane-wavefront  model  is  replaced  by  a  vector 
that  has  been  combined  over  all  the  modes,  we  might  know  the  individual 
modal  vectors  (which  depend  only  on  the  horizontal  wavenumbers  around  the 
array)  but  not  the  combining  weights. 

Hence  the  third  consequence  is  that  traditional  processing  methods  might  be 
ruled  out,  including  also  maximum-likelihood  and  maximum-entropy 
beamforming.  One  is  forced  to  such  methods  as  orthogonal  decomposition. 

The  fourth  consequence  is  that  the  combined  effect  of  a  larger  array  and 
refined  modelling  makes  the  processing  more  sensitive  to  non-stationari ty. 
As  an  example,  the  source  vector  of  the  normal-mode  model  (replacing  the 
steering  vector  of  the  plane-wavefront  model)  is  sensitive  to  changes  in 
geometry  caused  by  movement  of  source  and/or  receiver.  This  means  that  a 
non-stationary  cross-correlation  matrix  R  is  required,  making  it  difficult, 
perhaps  impossible,  to  obtain  a  good  estimate  of  this  central  quantity  in 
present  data-processi ng  methods. 
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Hence  the  fifth  consequence  is  that  processing  based  on  the  matrix  R  might 
be  abandoned  and  replaced  by  methods  that  take  into  account  the  modelling 
of  the  dynamic  nature  of  the  involved  processes.  One  obvious  candidate  is 
a  state/space  approach,  of  which  the  Kalman-Bucy  filter  is  the  best-known 
member. 

We  might  be  at  a  crucial  point  in  research  on  (passive)  sonar  data 
processing,  a  point  where  data-processi ng  researchers,  sound-propagation 
model  researchers,  and  research  managers  need  to  work  coherently  together 
for  the  common  aim  of  improving  passive  sonar  system  performance. 
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APPENDIX  A 

PARAMETER  ESTIMATION  TECHNIQUE  FOR  THE  ORTHOGONAL  DECOMPOSITION  METHOD 


Let  Rg,£=0,  1,  ...L,  be  constant  Hermitian  matrices: 

R  =  B„  -  lb*  =  \  AiEjE’  ,  L<N  (Eq-  A-  1 ) 

where  A,.  and  E.  ,  1=1,  ...  N,  are  eigenvalues  and  eigenvectors  of  the 

Hermitian  matrix  R  .  Therefore  all  the  eigenvalues  are  real  and,  by 
convention,  ordered  by  decreasing  value. 

The  problem  is  to  minimize  J  : 


N 

J  =  l  A2.  P  <  N  (Eq.  A. 2) 

i=N-P+l  1 


with  respect  to  the  parameter  vector  B: 

B  =  [b  ,  b  .  .  .  b.  ]T  .  (Eq.  A. 3) 

1  2  L 


One  iterative  algorithm  for  this  non-linear  problem  is: 
B(n+1)  =  B(n)  -  a(V|  J  +  MI)_1  vBJ 


(Eq.  A. 4) 


where  the  gradient  of  J  with  respect  to  B  is: 


G  =  V  * 


(Eq.  A. 5) 


and  the  so  called  Hess  matrix  of  J  with  respect  to  B  is: 


H  =  4bj  = 


a^j  a2J 

ab1ab1  • • • ab.ab. 

a2J  afj 

abLab1  ' ' ' abLabL 


(Eq.  A. 6) 


'  \ 
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Our  purpose  in  this  appendix  is  to  find  expressions  for  the  elements  of  G 
and  H.  We  make  the  following  definitions: 


D  (  )  =  1 

nr  '  9b 

m 


(Eq.  A. 7) 


p.  .  = E*D  (E.  )  , 

imj  i  mv  j 


(Eq.  A. 8) 


Y .  •  =  E*.R  E . 
imj  l  m  j 


(Eq.  A. 9) 


Using  the  operator  Eq.  A. 7  on  Eq.  A. 2  yields: 


°m(J)"2  1  *iVV  > 

i=N-P+l 1  1 


(Eq.  A. 10) 


Dk(0m(J))  =  2i=N!p+iDk(Xi)  W  +  Xi  ' 


(Eq.  A  11) 


We  now  use  many  of  the  properties  of  a  Hermitian  matrix  and  its  eigenvalues 
and  eigenvectors: 


(R  -  A. I)  E.  =  0 


Using  Eq.  A. 1  inEq.  A. 12,  and  operati ng  wi th  D  ,  gives: 

<•  WV  »  Ei  *<R- V>Dm<Ei>  =  0-  ’ 

Premul tiplyi ng  wi  th  E*j  gives: 


(Eq.  A. 12) 


(Eq.  A. 13) 


'EjRr»Ei  •Dm<Ai>EjEi4‘E|R-AiEj>»„(Ei)i0  • 


(Eq.  A. 14) 


A  Hermitian  matrix  has  orthonormal  eigenvectors  and  real  eigenvalues: 


-EjVrDi*1>  3ij  +  (Aj'  V  EjDm(Ei)  =  0  • 


(Eq.  A. 15) 


For  i  =  j: 


W  =  -  EflnEi  • 


(Eq.  A. 16) 


Dk<VAi»  =  -D^*)R»,El  -E’WEj)  ■ 


(Eq.  A. 17) 
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The  last  equation  (Eq.  A. 25)  has  to  be  manipulated  quite  a  lot: 


WJ))  2  i=N-P+1VikiYimi 


N  N-P 
+  2  I 
i=N-P+l 


2  \ . B  .  -  .y  .  .  '  X  -B ..  .v.  . 

j=1  i  ikj  jmi  r  jki *imj 


N  N 


2  I 
i=N-P+l 


2  A  .6 . ,  .Y  .  .  -  A  .B..  .y  . 
v  i kj ! jmi  j H l kj  *  j 


j=N-P+l 


jmi 


where  the  indices  i  and  j  have  been  exchanged  in  the  last 
last  line.  This  can  be  done  because  of  symmetry.  Using 
Eq.  A. 26  yields 


N 

D.  (D  (J))  =  2  2  y..  .y.  • 

K  lTl  hl/i  hmi 


i=N-P+l 


i  ki * i mi 


N  N-P  A. 

+  2  2  2  — (ViklYllB1+Ylkl-Yimi) 

■j=l^-p+^  j=l  j  i KJ  jroi 


N 

+  2  2 


i=N-P+l  j=N-P+l 


^Yikj  Yiki^ij^  Y.imi  ’ 


ij  jmi 


Dk(Dm(J)) 


WJ»  = 

where 


N  N-P  2  X. 

I 

i=N-P+l 


N-r  6  A. 

jll  "  Ai  Re*YikjYjmi* 


N  N 

•  1  1  2  YikiYimi  > 

i=N-P+l  j=N-P+l  1KJ  Jmi 


N 

2 

i=N-P+l 


N 


2  v .  .  Re?Y- 1  -Y  •  • }  . 

ij  1 ’l kj  *jmi ’  ’ 


4  A. 


v .  .  = 


ij  A.  -A,* 
=  2 


j  <  N-P 
j  >  N-P 


One  might  have  to  increase  the  eigenvalues  A.,  j=l,  ...  N-P, 
order  to  be  able  to  divide  by  Aj  *  A^ ,  j  <  N-P,'*  i  >  N-P. 


(Eq.  A. 26) 


term  of  the 
Eq.  A.  18  in 


(Eq.  A. 27) 


(Eq.  A. 28) 


(Eq.  A. 29) 


slightly,  in 


40 


SACLANTCEN  SR-84 


APPENDIX  B 

GENERALIZED  GONIOMETER  TECHNIQUE 


Let  R  be  a  Hermitian  matrix: 


R  -  F  F*  =  I  Fi  F*i  =  1  °f  °iAiA*D*  ’ 

i=l  i=l 


(Eq.  B.l) 


where  F.  ,  i=l,  ...  T  are  source  vectors: 
1  ’ 


F .  =  ct.  Z  a  .D  .=  o.D.A.  , 
1  l  t  m,  i  m,  l  ill’ 

m=l  ’  ’ 


(Eq.  B. 2) 


F  =  [Fi  ...  F, 


(Eq.  B. 3) 


D  .is  the  mode  vector  for  mode  m  and  source  i 

m,  i 


D  .  = 

m,  i 


’^kmxlCos(M" 


-  ik  x..cos(B . ) 
J  m  N  Vli ' 


(Eq.  B. 4) 


D.  =  [D,  .  .  .  .  D  . 
l  l,i  m,  l 


(Eq.  B. 5) 


and  A.  is  the  normalized  mode  amplitude  vector: 


Ai  =  Ial,i  * ‘ -  aM,i^ 


(Eq.  B. 6) 


A  i  A  i  =  1  ' 


(Eq.  B. 7) 


We  develop  some  properties  for  the  case: 


T  <  N 


(Eq.  B. 8) 


M  ^  N-T  . 


(Eq.  B.9) 


The  matrix  R  has  N  orthonormal  eigenvectors  En  ,  and  N  real  non¬ 
negative  eigenvalues  \n  ,  n=l  ...  N.  The  eigenvalues  are  indexed  in 
decreasing  order.  From  Eq.  B.7  it  follows  that: 


\  =  0 
n 


n=T+l,  ...  N 


(Eq.  B. 10) 
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such  that: 


R  EDADE0 


(Eq.  B. 11) 


where: 


E0  =  IE2  ...  Ey]  , 

'  \x  0  . . .  0  ■ 
0  X2  .  .  0 


(Eq.  B. 12) 


(Eq.  B. 13) 


[00...  xTJ 

Eq  spans  a  T-dimensional  subspace  Sq  .  The  source  vectors  ,  i=l,...T, 
will  be  linear  combinations  of  these  eigenvectors: 


F .  =  I  z  -  a  E  =  En  Z .  ,  1=1,  . . .  T 

1  r\=l  n»1  n  n  0  D  t  * 


(Eq.  B. 14) 


where: 


Zi  lzl,i  ZT ,  i 


(Eq.  B. 15) 


We  can  write  Eq.  B.14  more  compactly  as 


F  =  E0  V  z  . 


where: 


(Eq.  B. 16) 


ZMZ,...  ZT] 


(Eq.  B. 17) 


Inserting  Eq.  B.16  into  Eq.  B.l  gives: 


k  *  * 

R  =  ed  ad  2  Z  aJ  Ed  . 


(Eq.  B. 18) 


This  will  coincide  with  Eq.  B.ll  only  if: 


Z  Z*  =  I 


(Eq.  B. 19) 
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As  the  T  source  vectors  defining  F  are  independent  (discrete  sources), 
Z  is  non-singular,  and  has  an  inverse: 

Z'1  =  Z*  ,  (Eq.  B. 20) 

Z*  Z  =  I  .  (Eq.  B. 21) 

Hence  Z  has  to  be  an  unitary  matrix,  see  also  <6>. 

The  remaining  eigenvectors  will  define  an  (N-T)-dimensional  subspace  S 
that  is  orthogonal  to  the  source  vectors: 

E*  F.  =  0  ,  1=1,  .  .  .  T  (Eq.  B.22) 

where: 

Ex  =  [Em  ...  Em)  .  (Eq.  B.23) 

Substituting  Eq.  B.2  into  Eq.  B.22  yields: 

Ef  D.  A.  =  0  ,  1=1,  ...  T  .  (Eq.  B.24) 

Let  R  be  a  N-dimensional  hermitian  matrix  of  rank  T  .  We  want  to  match 

it  as  well  as  possible  to  a  matrix  R  given  by  the  model  of  Eq.  B.l.  The 

parameters  to  be  estimated  are  ,  and  A.  ,  i  =1 ,  ...  T. 

For  K  values  of  p  (assuming  p  is  properly  quantized),  we  try  to 
satisfy  Eq.  B.24,  using  eigenvectors  from  R.  Hence  we  want  to  minimize: 

Jk=|E*DkAk|2  k=l,  ...  K  (Eq.  8.25) 

with  respect  to  Ak  ,  while  at  the  same  time  obeying  Eq.  B.7.  The 
solution  is  the  eigenvector  of  D£  Et  E*  Dk  with  the  smallest  eigenvalue 

ek  : 

D£  lx  EJ  Dk.  Ak=  ek  Ak  ,  k=l,  ...K  (Eq.  B.26) 

Jk  =  ek  ,  k=l,  ...  K  .  (Eq.  B.  27) 

Let  us  plot  Jk  as  a  function  of  Pk  ,  and  assume  (and  hope)  that  this 

plot  has  T  minimal  points  (and  hope  that  they  are  close  to  zero).  These 

points  define  our  estimates  p^  and  A^  ,  1=1 ,  ...  T.  This  method  may  be 

called  a  generalized  goniometer  technique,  see  <7>.  One  way  to  obtain  an 
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estimate  ,  i=l,  ...  T,  is  to  combine  Eqs.  B.2  and  B.14: 

°1  °i  5i  -  Eo  kd  h  ■ 


_ L  _  * 

Premultiplying  Eq.  B. 28  with  Eg  gives 


2i  =  K  °i  *i  51 


From  Eq.  B.21  it  follows  that  should  have  unit  length: 

*  *  *  -1  _  * 

i  =  Z(  Z,=A(  0(EDA0E0  D,  A,  5^,  , 


[5*DD?DrJ  E*0^]'1  •  '*l—T 


Another  approach  is  the  following: 

T  T 

R  =  I  a?  D.  A.  A*.  D*.  =  I  a2.  R  .  . 

.  ,  11111  .,11 


(Eq.  B. 28) 

(Eq.  B. 29) 

(Eq.  B. 30) 

(Eq.  B. 31) 

(Eq.  B. 32) 


Then,  a.  ,  i=l,  ...  T,  can  be  estimated  by  minimizing  the  norm: 

J  =  |  R  -  R  A  =|  R  -  1  a2.  R.  |  (Eq.  B.33) 

i=l 

with  respect  to  a-  ,  i=l,  . . ,  T. 
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